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ABSTRACT 



Using the (3+1) formalism in general relativity, we perform the post-Newtonian(PN) approx- 
imation to clarify what sort of gauge condition is suitable for numerical analysis of coalescing 



X 

compact binary neutron stars and gravitational waves from them. We adopt a kind of transverse 



gauge condition to determine the shift vector. On the other hand, for determination of the time 
slice, we adopt three slice conditions(conformal slice, maximal slice and harmonic slice) and dis- 
cuss their properties. Using these conditions, the PN hydrodynamic equations are obtained up 
through the 2.5PN order including the quadrupole gravitational radiation reaction. In particu- 
lar, we describe methods to solve the 2PN tensor potential which arises from the spatial 3-metric. 
It is found that the conformal slice seems appropriate for analysis of gravitational waves in the 
wave zone and the maximal slice will be useful for describing the equilibrium configurations. The 
PN approximation in the (3+1) formalism will be also useful to perform numerical simulations 
using various slice conditions and, as a result, to provide an initial data for the final merging 
phase of coalescing binary neutron stars which can be treated only by fully general relativistic 
simulations. 



1. Introduction 



Kilometer-size interferometric gravitational wave detectors, such as LIGO and VIRGO , 
are now in construction aiming at direct detection of gravitational waves from relativistic astro- 
physical objects. Coalescing binary neutron stars are the most promising sources of gravitational 
waves for such detectors. The reasons are that (1) we expect to detect the signal of coalescence 
of binary neutron stars about several times per year' 31 , and (2) the wave form from coalescing 
binaries can be predicted with a high accuracy compared with other sources' 11 . 

In the case when the orbital separation of each star is large compared with the radius of 
neutron stars, i.e., in the so-called inspiraling phase, binary neutron stars are evolving in the 
adiabatic manner due to gravitational radiation reaction with much longer time scale than the 
orbital period. As for the inspiraling phase, the theoretical investigation is usually done by 
the point particle approach using the PN approximation in general relativity [4 ' 5 ' 6,71 . Since the 
separation is large compared with the neutron star radius, the hydrodynamic effect is small 
enough and we can regard each star of binary as a point particle. Theoretical studies for 
such a phase is potentially important because by comparing the observational signal with the 
theoretical prediction of the signal of inspiraling binary, we may be able to know not only the 
various parameters of binary [8,91 , but also the cosmological parameters' 101 . 

After a long time emission of gravitational waves, the orbital separation becomes comparable 
to the radius of the neutron star. Then, each star of binary neutron stars begins to behave as 
a hydrodynamic object, not as a point particle, because they are tidally coupled each other. 
Recently, Lai, Rasio and Shapiro' 111 have pointed out that such a tidal coupling of binary neutron 
stars is very important for their evolution in the final merging phase because the tidal effect causes 
the instability to the circular motion of them. Also important is the general relativistic gravity 
because in such a phase, the orbital separation is larger than ~ 10% of the Schwarzschild radius 
of the system. Thus, we need not only a hydrodynamic treatment, but also general relativistic 
one to study the final phase of binary neutron stars. 

Fully general relativistic simulation is sure to be the best method, but it is also one of the 
most difficult ones. Although much effort has been focused and much progress can be expected 
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of the main reasons is that we do not know the behavior of the geometric variables in the 
strong field around coalescing binary neutron stars. Owing to this, we do not know what sort 
of gauge condition is useful and how to give an appropriate general relativistic initial condition 
for coalescing binary neutron stars. 

The other reason is a technical one: In numerical relativistic simulations, gravitational waves 
are generated, and in the case of coalescing binary neutron stars, the wavelength is the order 
of A ~ vr J R 3 / 2 M" 1 /2 ) wnere R anc i M 

are the orbital radius and the total mass of binary, 
respectively. Thus, we need to cover a region L > A oc R?l 2 with numerical grids in order 
to perform accurate simulations. This is in contrast with the case of Newtonian and/or PN 
simulations, in which we only need to cover a region A > L > R. Since the circular orbit of 
binary neutron stars becomes unstable at R ^ 10M owing to the tidal effects lnl or the strong 

[121 

general relativistic gravity , we must set an initial condition of binary at R ^ 10M. For such a 
case, to perform an accurate simulation, the grid must cover a region L > A ~ 100M in numerical 
relativistic simulations. When we assume to cover each neutron star of its radius ~ 5M with 
~ 30 homogeneous grid' 14 ' 1 '' 1 , we need to take grids of at least ~ 500 3 , but it seems impossible 
to take such a large amount of mesh points for the present power of supercomputer. At present, 
we had better search other methods to prepare an initial condition for binary neutron stars. 

In the case of PN simulations, the situation is completely different because we do not have 
to treat gravitational waves explicitly in numerical simulations, and as the result, only need to 
cover a region at most L ~ 20 — 30M. In this case, it seems that ~ 200 3 grid numbers are 
enough. Furthermore, we can take into account general relativistic effects with a good accuracy: 
In the case of coalescing binary neutron stars, the error will be at most ~ Mj R ~ a few x 10% 
for the first PN approximation, and ~ (M/R) 2 ~ several % for 2PN approximation. Hence, if 
we could take into account up through 2PN terms, we would be able to give a highly accurate 
initial condition(the error ^ several %). For these reasons, we consider the 2.5PN hydrodynamic 
equations including 2.5PN radiation reaction potential in this paper. 

The purpose of this paper is twofold: One is to establish the basic formulation of the 2.5PN 
hydrodynamic equation, and the other is to investigate what kind of gauge condition is appro- 
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from them. As for the PN hydrodynamic equation, Blanchet et al. have already obtained the 
(1+2.5)PN formula' 161 . Although their formula was very useful for PN hydrodynamic simula- 
tions including the radiation reaction [14,lj ' 17] ; they did not take into account 2PN terms. Also, in 
their formula, they fixed the gauge conditions to the ADM gauge, but in numerical relativity, it 
has not been known yet what sort of gauge condition is suitable for simulation of the coalescing 
binary neutron stars and estimation of gravitational waves from them. For these reasons, we 
shall investigate several gauge conditions using the (3+1) formalism in general relativity. 

This paper is organized as follows. In section 2 we present the (3+1) formalism of the 
Einstein equation and the equations for the PN approximation. Several slice conditions are 
imposed in section 3. In section 4, the quadrupole radiation-reaction potential is calculated in 

[18] 

combination of the conformal slice and the transverse gauge. It is found that this combination 
of the gauge conditions simplifies the calculation of the back reaction potential. The methods to 
solve the 2PN tensor potential are discussed in detail for the sake of actual numerical simulations 
in section 5. The conserved energy and conserved linear momentum are refered in section 6. 
Section 7 is devoted to summary. 

We use the units of c = G = 1 in this paper. Greek and Latin indices take 0,1,2,3 and 
1, 2, 3, respectively. 

2. (3+1) Formalism for Post-Newtonian Approximation 

2.1. (3+1) Formalism 

We consider the (3+1) formalism to perform the PN approximation. In the (3+1) formalism, 
the metric is split as 



(2.1) 



and 



n M = (-a, 0), 




(2.2) 
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respectively. Then the line element is written as 

ds 2 = -(a 2 - i3ift)dt 2 + 2(3 i dtdx i + lij dx i dx j . (2.3) 

Using the (3+1) formalism, the Einstein equation 

GV = SvrTV, (2.4) 

is split into the constraint equations and the evolution equations. The formers are the so-called 
Hamiltonian and momentum constraints which respectively become 

tiR - K i:j K ij + K 2 = Wttph, (2.5) 

DiK { - - DjK = 8tcJj, (2.6) 

where K^, K, trR and D{ are the extrinsic curvature, the trace part of K^, the scalar curvature 
of 3D hypersurface and the covariant derivative with respect of 7^. pu and Jj are defined as 

PH = T^n", 

(2.7) 

Jj = -t^Yj- 

Evolution equations for the spatial metric and extrinsic curvature are respectively 

J^ 7 y = -2aKij + Dify + Dj(3i, (2.8) 

= a(R i:j + KKtj - 2KuK l A - D.Dja 
M _ (2.9) 

+ (Djf3 m )K m i + (Di(3 m )K mj + f3 m D m K ij - Ma[s i3 + - 7 y(p# - S 1 ^, 
^ 7 = 2 7 (-«^ + A/3 i ), (2.10) 

^-K = a(trR + K 2 ) - DWia + ft DjK + Ana(S\ - 3p H ), (2.11) 

where Rij, 7 and Sij are, respectively, the Ricci tensor with respect of 7^, determinant of 7^ 
and 

^ = TH 7 W r (2.12) 



To distinguish the wave part from the non-wave part (for example, Newtonian potential) in 
the metric, we use 7^ = ip Hj instead of 7^. Then det(7y) = 1 is satisfied. We also define Aij 

as 

Aij = V~ 4 Ay ee ^ 4 (iYy - \lijK). (2.13) 

We should note that in our notation, indices of are raised and lowered by 7^, so that the 
relations, A* ■ = A* ■ and A 1 - 7 = ■0 4 A y , hold. Using these variables, the evolution equations 
(2.8)-(2.11) can be rewritten as follows; 

d _ n j- _ 9/3' _ 2_ 9/3' . 

- 7ij = -2o^ + 7^ + 7,7^- - 37^^, (2.14) 



9 ~ 1 
dn~ Aij -ft ° 



<9/3 m ~ <9/3 m ~ 2 <9/3 m ~ a / 1 

+a(KA ij - 2A a A l ? ) + -£-A mj + -£-A mi - ~^~Mi - ^{Sij ~ 3^, 



(2.15) 

d ^ ( K ®@ % \ (2 iq) 
dn 6 V dx % J ' 

|-K = + - ^Aa - 47*ty,*a i + 47ra(S\ + (2.17) 
where .D, and A are the covariant derivative and Laplacian with respect to 7^ and 

9 d ■ d 

dn dt ^ dx l ^ ^ 
The constraint equations are also written as 

M = - 2n PH ^ 5 - y { A » Aij - I r2 ) » ( 2 - 19 ) 



and 



DM &A \) - l^ 6 D t K = Sn^Ji, (2.20) 



whprp trR is thp sralar mirvfltiirp with rpsnprt tn 



Now let us consider Rij in Eq.(2.15), which is one of the main source terms of the evolution 
equation for A{j. First we split Rij into two parts as 



Rij — Rij + R-fj, 



(2.21) 



where Rij is the Ricci tensor with respect to 7^, 



R% = ~DiDj^ - ^jD k D k ^ + ^{D^){b^) - ^ij(D k i;)(D k ^. 



(2.22) 



Using the property of det^j) = 1, Rij is written as 



R 



l H {llj,ik + llijk ~ lij,lk) + l K \k(llj,i + 1H,3 ~ 



~ 1 kj L li 



(2.23) 



where ^ denotes d/dx l and f^- is the Christoffel symbol with respect to %j. We split 7^ and 



7 y as 5ij + hij and <5 y + / y , where 5ij denotes the flat metric, and rewrite Rij as 



R 



-Kj,kk + hji : ii + hu ij + f k \ k (hij :i + h U j - hijj) 



+ f kl {h k j^i + hkiji - hij yk i) 



~ 1 kj L li 



(2.24) 



In this paper, we consider only the linear order in hij and assuming \hij\, 1. (So 

that, hij = — P 3 .) Such an assumption is justified because in this paper, we choose a gauge 
condition, in which hij is a 2PN quantity(see below). This implies that we neglect higher PN 
effects such as the non-linear coupling between gravitational waves themselves, but does not 
imply that we neglect the non-linear coupling between the Newtonian potentials themselves and 

between gravitational waves and the Newtonian potentials. In other words, although we can not 

[21] 

see the non-linear memory of gravitational waves , we can see the tail term of gravitational 
waves and can derive the exact quadrupole formula(see below). Here, to guarantee the wave 
property of 7^, we impose a kind of transverse gauge* to hij as 



(2.25) 



This condition is guaranteed by (3 l which satisfies 



-^j^k = (-2aA tj + inP'j + ~ hijP l ,i) • (2-26) 

Using the above conditions, Eq. (2.24) becomes 

i^ = ~A /Iot /^- + C>(/i 2 ), (2.27) 

where A fiat is the Laplacian with respect to Sij. Note that tri? = 0(h 2 ) is guaranteed in the 
transverse gauge because the traceless property of hij holds in the linear order. 

Finally, we show the equations for the perfect fluid. The energy momentum tensor for the 
perfect fluid is written as 

T» v = (p + p£ + P)u^u v + PgW, (2.28) 

where w M , p, e and P are the four velocity, the mass density, the specific internal energy and the 
pressure. The mass density obeys the continuity equation 

V>^) = 0, (2.29) 

where V M is the covariant derivative with respect to g^ v . The explicit form is 

where p* is the conserved density defined as 

p* = aip 6 pu°. (2.31) 
The equations of motion and the energy equation are derived from 

V M T^ = 0. (2.32) 

Explicit forms of them become 

dSi , d(s^) = _ - _ + ■ _ i jk (2 ^ 33) 



and 



where 



5< = c# 6 (p + + P)u°ui = p*(l + e + ^)«»(= ^ 6 ^i). 

g o = a0 8 (p + ^ + P)(« O ) 2 (= (Pg+ / ) ^ ), 
if = aip 6 peu° = p*e, 



(2.35) 



Finally, we note that in the above equations, only (5 l appears, and does not, so that, in the 
subsequent section, we only consider the PN expansion of not of 

2.2. Post-Newtonian approximation 

Next, we consider the PN approximation of the above equations. First of all, we review the 
PN expansion of the variables. Each metric variable, which is relevant to the present paper, is 
expanded as 

i> = 1 + (2)^ + (4)^ + (6)^ + (7)^ + • • • , 
a = 1 + [2) a + (4) a + (6) a + (7) a + • • • ) 

= l-U+i K — + Xj+ (6) a + (7) a + . . . , 
P = (3)A + (5)A + (6)A + (7)A + (8)A + • • • , (2.36) 
= (4) hij + (5) hij + ... , 

= (3) A? + (5) A? + (6)Aj + • • • , 
# = (3)^ + (5)# + (6)# + •••, 

where subscripts denote the PN order(c~ n ) and U is the Newtonian potential satisfying 

A flat U = -4np. (2.37) 



X denends on the slice condition, and in the standard PN raiiffe . we nsnallv use <J> = 



-XII 



which satisfies 



Aflat® 



rr 1 3P \ 



(2.38) 



Note that the terms relevant to the radiation reaction appear in ^ip, {j)Ct, (8)A and (5)/iy, and 
the quadrupole formula is derived from and (5)/iy. 

The four velocity is expanded as 



M = 



^ + ^ 2 + (3)/V-x)+0(c- 6 ), 
1 + (\v 2 -U) + (jj, 4 + ^ 2 f/ + ^ + X)] + 0(c~ 6 ), 

=«' [i + (^ 2 + ^) + (jj« 4 + ^ + \u 2 + (3)/V - x)] + 0(c 

(5) A + (3)A(t^ 2 + 3C/) + (4)/lj. 



Ui =v" + 



+ 



+ U V2 

'3 .4 i 7 2, 



(2.39) 



where f 2 = v l v l . The PN expansion of the relation u^u^ = — 1 becomes 

(cm ) 2 = 1 + 7« WiUj - 

= 1 + v 2 + + 4v 2 f7 + 2 (3) A^ + 0(c~ 6 ). 
Thus Jj and SV/ are respectively expanded as 

1 + (v 2 + e) + j> + t; 2 (417 + £ + -) + 2 (3) /W} + 0( C - 6 ) 



(2.40) 



PH =P 



( V V + ^) + { (t; 2 + 6C/ + e + + ^ (3 )& + ^'(3) A + 2 ^-^} (2 ' 41) 

+ 0(c- 6 )", 

t; 2 + — + {2(3)^ + v 2 (v 2 + 417 + £ + -)} + 0(cT 6 ) 



5, =P 



■0(and a in the conformal slice) is determined by the Hamiltonian constraint. At the lowest 
order, it becomes 



A/Jot (2)^ = ~ 27r P- 



(2.42) 



Thus. to\a = —2<r,\ib = —U is satisfied in this DaDer. At the 2PN and 3PN orders, the Hamilto- 



nian constraint equation becomes, respectively, 

A/iat(4)^ = -2np{v 2 + e+ b -U), (2.43) 



and 



A/iat (6 )^ = - 27rp{t; 4 + t; 2 (e + - + — if) + 2 (3) /V + -eU + -t/ 2 + 5 (4) ^} 

+ ^(4) h ij U ,ij ~ g((3)Aj(3)^j - 3(3)^ 2 )- 



(2.44) 



The term relevant to the radiation reaction first appears in tj\i]) and the equation for it becomes 

^flat{7)^ = ^(5) h ijU,ij. (2.45) 

Hence, may be also relevant to the radiation reaction and whether it may or not depends 
on the slice condition. 

From Eq.(2.26), the relation between (3) Ay and ^(3i becomes 

2 

- 2 (3)Aj + ( 3 )&j + (3)/^,i - 3 <J v(3)A,i = 0. (2.46) 

(3) Ay must also satisfy the momentum constraint. Since (3) Ay does not contain the transverse- 
traceless(TT) part and only contains the longitudinal part, it can be written as 

(3)Ay = {3)WiJ + {3) W jti - -S ij{3) W Kk , (2.47) 

where (3)Wj is a vector on the 3D hypersurface and satisfies the momentum constraint at the 
first PN order as follows; 

A/ Io t(3)Wi + ~ {3) Wjji - - {3) Ki = 8irpv\ (2.48) 



From Eq.(2.46), the relation, 



i^3i = 2,^W; 



holds and at the first PN order, Eq.(2.16) becomes 



3U = - {s) K + {3) Pu , 



(2.50) 



where U denotes the derivative of U with respect to time, so that Eq.(2.48) is rewritten as 



&flat(3)Pi = ISTrpu 4 + ( (3) ^ - Uij. 



(2.51) 



This is the equation for the vector potential at the first PN order. 



From the next order, r n -\f3i is determined by the gauge condition, hijj = 0. Making use of 
the momentum constraint and the 2PN order of Eq.(2.16), 



6(4)^ - 3 (3) /W,/ - 7}U(2 i3) K + 3U) + {5) K = (5) A,I , 



(2.52) 



the equation for is written as 
Aflat(5)Pi = !6vrp 



v i (y2 + 2U + e + ^j+ (3) A] - 8f/, i(3 )iy 



+ (5)K,i - U^Kj + \u m K - 2(4)^ + \{UU\i + ( (3 )/TO,i • 



(2.53) 



Since Jj at the 1.5PN order vanishes, the merely geometrical equation for ( 6 )/3j is given by 



Aflat(6)f3i = (6) K ,i ■ 



(2.54) 



Then, let us consider the wave equation for h^. From Eqs.(2.14), (2.15), (2.21) and (2.27), 
the wave equation for hij is written as 

f d 2 d 2 \, 

□ hij ={l~ ¥ ) Aflathij + -Qp) h V 

2a r ~ (bibj - \zfijk)i> + ^ (b^bji, - l -^jb k i,b k ^) 



- [bibj - -7yAja - -{b^bja + bjipbia - -^jb k ^b k aj 



1> 

2n 2 [ KAj - 2A u A l j) + 2a[f3 m i A mj + (3 m jA mi - ^/3 m m iy) 

2 



<9n 



(9a 



(2.55) 



where 



2 



□ =- W + ^flat. (2.56) 

We should note that (4)7y has the TT property, i.e., (A) T ij,j — an d (4) T u — 0- This is a natural 
consequence of the transverse gauge, hijj = and ha = 0{h 2 ). Thus (4)/iy is determined* from 

Aflat(4)hij = (4) T ij- (2-57) 

(5)/iy is derived from 

(5)hij(t) = ^ y ( 4)Ty(t,y)rf 3 y, (2.58) 
and the quadrupole mode of gravitational waves in the wave zone is written as 



'■'* = -!■ hm /^ (t ,- |X 7 yl ' y) A. (2.59) 



In the subsequent section, we derive the quadrupole radiation-reaction metric in the near zone 
using Eq.(2.58). 

Finally, we show the evolution equation for K. Since we adopt slice conditions which do 
not satisfy K = 0(i.e., maximal slice condition), the evolution equation for K is necessary. The 
evolution equations appear at the 1PN, 2PN and 2.5PN orders which become respectively 

8 . , . ... /> 



dt 



(3) J 



Airp^v 2 + e + 2U + 3-) - A flat X, (2.60) 



— {5) K =4tt P \2v a + v 2 (6U + 2e + 2-)-[e + -)U- AU 2 + 4 (4) ^ + X + 4 (3) /^ 

+ (3)A?(3)A? + 3(3)^ _ {4)hijU,ij + ( 3 )A(3)-ftr,i (2-61) 

- ^UU k U k - U, k X :k + 2U m i>, k - A flat{6) a + 2UA flat X, 

d 

fa(6) K = ~^flat(7) a ~ (5) h ijU,ij- ( 2 -62) 

We note that for the PN equations of motion up to the 2.5PN order, we need ( 2 )Q;, (4) a, 
(6) a, (7) a, (2) ^, (4) ^, (3) A, (5) #, (6) /3j, (4) /iy, (5) /iy, (3) K, (5) K and (6) K. Therefore, if we solve 



* Since 0(h 2 ) turns out to be 0(c 8 ), it is enough to consider only the linear order of hij in the case when 



the above set of the equations, we can obtain the 2.5 PN equations of motion. Up to the 2.5PN 
order, the hydrodynamic equations become 



dSi d(SiV 



dt 



+ p* [U,i{l + e + ^ + h 2 -U+^ + Av 2 U + {^v 2 -u)(s + ?-)+ 3 (3) /^} 

( P V 2 \ 2 

-X ji \l+e + — + —)+2v (4) ^,i - ( 6 )Q!,i - (7)a i 
P v 2 



dH d(Hvi) 
dt dxi 



+ ^*((4)fyfc,t + (5) h jk ,i) + 0(c- 8 ) 

) + ^{(r 2+3C/ M +0 <^>] 



(2.63) 
(2.64) 



where we make use of relations 



aS° = p Jl + e + - + V — + V — (e + -) + \v A + 2v 2 f/ + + +O(o 



^ (l + e + + y + 3C/) + (3) A + 0(c- 5 ) 



(2.65) 



3. Slice Conditions 

[18] 

In this section, we perform the PN analysis using the conformal slice , maximal slice and 
harmonic slice' 201 which are often used in 3D numerical relativity. Among them, we find that the 
conformal slice seems most tractable and useful to estimate gravitational waves in the far zone, 
while the maximal slice is suitable for describing the equilibrium configurations. Hence, first 
of all we describe the property of the conformal slice and then mention the properties of other 
slices. 



3.1. CONFORMAL SLICE 



The conformal slice is defined 1181 as 



( n 2 3 2 <\ 

a = expl— 2e — -e — -e I, 
V 3 5/ 



(3.1) 



where e = ip — 1. In the conformal slice, becomes 



(2)« = - 2 (2)VS 

(4)« = 2 ((2)^) 2 - 2 (4)VS 

(6) « = - 2 ((2)^) 3 + 4(2)^(4)^ - 2 (6)^ ? 

(7) " = - 2 (7)^S 



(3.2) 



Although in the usual PN approximation we need to solve the Poisson equation for the lapse 
function, this slicing saves solving it. 

In the conformal slice, equations (2.15) and (2.17) are rewritten as 



-At 



1 a 



2a 



( 



DjDjt ■ - ^D k D k ^)Ul + e + e 2 + e 3 + e 4 ) 



(3 + 6e + 6e 2 + 6e 3 + 6e 4 + 12e 5 + 10e 6 + 8e 7 + 6e 8 + 4e 9 + 2e 10 )l (3.3) 



+ aiKAij - 2A il A l j ) + P^Amj + {5™~A mi - -(3 m m A tj 



l 

3 



and 



dK a 



2 ~ 



A^(l + t 1 + e 4 ) - -^D k ^D^(3e° + 2e b + 2e' + e 8 + e 9 ) 



(3.4) 



whprp wp nsp thp TT nrnnprtv as wpll as thp linpar flnnrrorirnfltinri fnr Yi.aa in thp ahnvp pnnatinn 



Then Eq.(2.55) is written as 



D ^=-(^- 1 ) A ^ + (^ ^ 



4a 

+ 



DiDjip - + e + 6 2 + 6 3 + 6 4 ) 

(3 + 6e + 6e 2 + 6e 3 + 6e 4 + 12e 5 + 10e 6 + 8e 7 + 6e 8 + 4e 9 + 2e 10 )] 
+ 2a 2 (KA tj - 2A ll A l j ) + 2a(f3 m i A mj + ^A mi - ^/3 m m iy) 



zT ij ) 



=(4)Ty+C(c b ), 

where we use e = 2£7 which holds in the PN approximation. 
In the conformal slice, the asymptotic form of e becomes 



(3.5) 



where we use e — tp — 1 and ip satisfies 

= -2ixp H ^ - ^ (AijA* - ^ 2 ) . (3.6) 

Eq.(3.5) is expanded as follows; 

□ hij =(uUij - ^SijUAfteU - 3UiUj + 8ijU k U k } 

- 167r(puV - U ljP v 2 ^ - fafcj + (3) 4i - pi j{3 )Pk,k) + 0(c- 6 ) ( 3 - 7 ) 

6x 



M ADM /q ox 

e ~ . (3.8) 



Thus a behaves as 1 — Madm/t at spatial infinity. This means that in the conformal slice, the 
metric at spatial infinity becomes the static Schwarzschild's one. This property seems helpful 
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Also, we have an advantage to derive the radiation reaction potential in this slice; From a 
relation ^)Ot = —2^ip and Eq.(2.45), we have 



A/jot (7)a = -(5) h ij(t)U,ij- (3.9) 



Thus, the radiation reaction potential is derived as 



(7)« = 4?r / Uijd a x. (3.10) 



Finally, we comment on the following week point of the conformal slice; in the conformal 
slice, the evolution equation for becomes 



{3) k = 4irp(v 2 + 3---u). (3.11) 



Since K does not vanish, K continues to change even in the case of a stationary spacetime. 
Thus, it seems inconvenient to describe equilibrium configurations of stars and binary systems 
in the conformal slice. To describe equiriblium configurations, we had better use the slice, such 
as the maximal slice, where K = is satisfied. 

3.2. Maximal Slice 

The maximal slice is given by 

K = 0, (3.12) 
and this equation leads to the equation for a as 

D k D k a = a(Aij£ j + 4n(E + S 1 ^. (3.13) 
At the first PN order, the equation becomes 

/ 3-P \ 

A f i a t {4) X M s = 4vrp [2v 2 + e + — + 2U), (3. 14) 
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relation holds; 

X cs = -2 {4) ^, (3.15) 
where CS similarly denotes "conformal slice". Using the above equation, we rewrite Xms as 



X M s = -2 ( 4)^ + Y. (3.16) 

Then the equation for Y becomes 

A f i at Y = 4n(pv 2 + 3P- l - P u). (3.17) 



[23] 

We should also note that by means of the virial theorem , the integration of the source 



term for Y can be written as 



where 



■ (t) = J px i x j d 3 x. (3.19) 
Hence, the behavior of Y far from the matter becomes 



2r 

In total, the behavior of a in the wave zone becomes 



Y ~ -^ki- (3.20) 



a ~ 1 



\{M + \h\ (3.21) 



Therefore, contrary to the conformal slice, in the maximal slice, the spurious time-dependent 
term is included in a in the wave zone. Since the metric does not approach the static Schwarzschild| 
metric even in spatial infinity, the maximal slice is inconvenient to distinguish a wave part from 

n i~\ ft tx ro tra no "vf o on r»h o o "f" h a l\Toii7+" nn ion nAranf i o 1 



In the case of the maximal slice, the equations for the shift vector are obtained by simply 
taking K = in Eqs.(2.51) and (2.53). Also, it is found that the equation for is the same 
as that in the conformal slice: The right-hand side of Eq.(3.13) has no 0(c~ 7 ) terms. Therefore 
Eq.(3.13) becomes 

Aflat(7)® = -(5) h ij U ,ij- ( 3 - 22 ) 

Finally, we show the wave equation for hij in the maximal slice as 
□ hij = - 2(1%- - ^A /Iot y) + (uUij - ^SijUA flat U - 3UjU d + 6ijU k U k ) 
- 167r(puV - ^S ijP v 2 ^ - (( 3 )/3ij + ( 3 )% - g<%( 3 )/\fc) + 0(c- 6 ). 



(3.23) 



3.3. Harmonic Slice 

The condition for the harmonic slice is 

□ t = 0, (3.24) 

which becomes in the (3+1) terminology, 

a + a 2 K - p { a,i = 0. (3.25) 

Differentiating this equation with respect to time, the wave equation for the lapse function is 
derived as 

□ a =^a\s\ - 3p H ) ~ A - A flat )a - ^-&^D ia - ^A^ 



+ 2aaK + —R + a 3 K 2 + a 2 (?DiK - fra^ - (3 i a i 



where A Q is expanded as follows, 

,jj2 



A a = 4tt P [i + (y 2 + 3- - -)] + A /Iot (— - 2 (4) V) + 0(c" 6 ). (3.27) 



P 

This equation is formally solved by using the retarded Green function and the Taylor expansion. 
For example, we obtain the Newtonian and first PN order lapse function 



(2) 



a 



= - f**£=s = -\{'** +0 <r-*), (3-28) 



( 4)« = -\J d s yp\* - y| - / d 3 y (^!+^_jpg + l ^ _ w . (3 29) 

Thus, at the spatial infinity, we find the following behavior 

(4)0! + 2 (4) ^ ~ (/fcfc - ^nV/fci) , (3.30) 

where = x l jr. From these equations we find that at the spatial infinity the lapse function does 
not behave as 1 — M/r + 0(r~ 2 ) unlike in the conformal slice, but behaves as 1 — (M + T(t))/r + 
0(r~ 2 ). Thus the harmonic slice is also inconvenient to distinguish a wave part from non-wave 
parts. 

The quadrupole radiation reaction potential takes the following rather lengthy form. 

i r h TT (3.31) 

1 / ,Q (^hijUi 



47TJ x-y 



This expression is similar to Chandrasekhar's one in the harmonic gauge 1241 and indicates that 
the fifth time derivative of the quadrupole moment appears in the reaction force, which is not 
convenient to treat in numerical calculations. 

4. The Radiation Reaction due to Quadrupole Radiation 

This topic has been already investigated by using some gauge conditions in previous pa- 
pers' 24 ' 19,161 . However, if we use the combination of the conformal slice and the transverse gauge, 
calculations are simplified. This is why we briefly mention the derivation of the radiation reaction 
potential in this section. 

In combination of the conformal slice and the transverse gauge, Eq.(2.58) becomes 

(5)M*) J[-l^{pV l V 3 - T^yV) 

+ (iJUij - U tj UA flat U - WiUj + SijUkUk)] d 3 y (4.1) 
1 d f / ■ ■ 2 • \ 3 

+ ~4^Qt J [w^j + ( 3 )^ ~ 2 S i3wPk,it) d y- 

From a straightforward calculation, we find that the sum of the first and second lines becomes 
and the third line becomes 6^-^/5, where 4-^) = d 3 4-a/dt 3 . (This calculation is 



replaced by a fairly simple one noticing the transverse property of (4)Ty. It is described in the 
appendix A.) Thus, (5)/%- in the near zone becomes 

(5)hij = (4.2) 

where 

■irij = hj - -kjhk- (4.3) 

Since hij has the transverse and traceless property it is likely that (5)^ remains the same for 
other slices. However it is not clear whether the TT property of is satisfied even after the 
PN expansion is taken in the near zone and, as a result, whether t^hij is independent of slicing 
conditions or not. The fact that slicing conditions never affect (tyhij is understood on the ground 
that (4)Ty does not depend on slices, which will be shown in the section 5. 

Then the Hamiltonian constraint at the 2.5PN order turns out to be 

^fiatii)^ = -U$U tij = A flatX ,ij , (4.4) 



where \ is the superpotential [221 and defined as 



X = -J p|x-y|d 3 y, (4.5) 
which satisfies the relation AfiatX — —2U. From this, we find (7\i[) takes the following form, 

1 av f „.( xJ -y j )^ 

(4.6) 



1 jlW(_„3tt. , f PW 



,3 



Therefore, the lapse function at the 2.5PN order, = —2^ip, is derived from U and U r , 
where U r satisfies' 161 

AUr = _ M f) p . x 3\ (4.7) 

Since the right-hand side of Eq.(2.62) cancels out, (6)-^ disappears if the does not exist on 
the initial hypersurface, which seems reasonable under the condition that there are no initial 
gravitational waves. Also, ( 6 ) A vanishes according to Eq.(2. 54). Hence, the quadrupole radiation 
reaction metric has the same form as that derived in the case of the maximal slice' 19 ' 16 '. 



From Eq.(2.33), the PN equations of motion becomes 



v { + M j = ~ + U ,i + F l PN + F ? PN + F ?- 5PN + °( c ~ 8 )> ( 4 - 



where Fl PN and Ff PN are, respectively, the 1PN and 2PN forces and conservative ones. Since 



the radiation reaction potentials, f^hij and moi, are the same as those by Schafer 1191 and Blanchet, 
Damour and Schafer [161 in which they use the ADM gauge, the radiation reaction force per unit 

?2.5PN 



mass, F^^ N = F[, is the same as their force and 

F i = -{i(5)hijv j y + v k v J M5) h ij + {7) a^ 

T 4 ^r(3) jy, 4 _e(3) k j 2 ( 3 ) d f (x k -y k )(x l -V) ^ 3 

-(^,-V) + v k + -F ( '— / p(t,y) d y 



Since the work done by the force (4.9) is given by 



(4.9) 



W = J pF[v i d 3 x 
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(4.10) 



we obtain the so-called quadrupole formula of the energy loss by averaging Eq.(4.10) with respect 
to time as 



5. Strategy to obtain 2PN tensor potential 

In this section, we describe methods to solve the equation for the 2PN tensor potential u\hij. 
Although Eq.(2.57) is formally solved as 

u)hij(t,-x) = -— ^ r d s y, (5.1) 

4vr y |x — y| 

it seems difficult to estimate this integral in practice since u\Tij — > 0(r -3 ) for r — > oo and 
the integral is taken all over the space. Thus it is desirable to replace this equation by some 
tractable forms in numerical evaluation. In the following, we show two approaches: In section 
5.1, we change Eq.(5.1) into the form in which the integration is performed only over the matter 
distribution like as in the Newtonian potential. In section 5.2, we propose a method to solve 

F,n (0. n« the> Vinnnrlnvv vnlnp nrnhlpm 



5.1. Direct integration method 



The explicit form of (4)Ty is 

(4) ry = - 2dij(X + 2 (4) ^) + C/^C 7 - 3C Wj + SijU k U k - 167r(puV - ^-pu 2 ) 

(3)Aj + (3)/fy,i - 2^i(3)/5fc,fcJ, 



where 

d 2 1 

Although u\Tij looks as if it depends on the slice condition, the independence is shown as follows. 
Eq.(2.51) is rewritten as 

Aflat(3)Pi = A flatPi + (3) K ,i, (5-3) 

where 

Pi = -4 / ^-fy -\{Jp\*- y\ d3 y) t ■ ^ 



This is solved as 



^"-M/js^X- (5 ' 5) 



From Eqs.(2.43) and (2.60), we obtain 



(3)^ = -A /Iot (X + 2 (4) ^) + 4np(y 2 + 3^ - |). (5.6) 



Combining Eq.(5.5) with Eq.(5.6), the equation for ( 3 )/3j is written as 



U2 + 3P _ pU/2 \ 

(3) A = « - (X + 2 (4) ^),i - { / ^ -— ^ 3 y| . . (5.7) 



Using this relation, the source term, (4)7y, is split into five parts 



(S) , (£/) , (C) . (p) . 00 /r Q x 



where we introduce 



(4)^ =UU,ij - ^5 tj UA flat U - SUiUj + 5ijU k U k , 

«> T « = 4 dJ J W^J\ d y + ^ J W^7\ d y - 3** ft? J W^7\ d y ' (5.9) 



if =dij J p|x-y|A, 



(4) r ii 
(4) 



Thus it becomes clear that (4)/iy and (5)/iy as well as (4)Ty are expressed in terms of matter 
variables only and thus do not depend on slicing conditions. 

Then, we define A //at(4) /i^ } = {a) t^\ A/^/ig = ^r^, A /Zat(4) /^ C) = (4) ^ C) , 
^flat(^)h\f = (4) T if and A// at ( 4 )/i^ = (±)t\J\ and consider each term separately. First, 
since {A) T ij ls a compact source, we immediately obtain 

(4)^f = 4 / V |x _ y| L *v- 

Second, we consider the following equation 

A /Iot G(x,yi,y 2 ) = , p r ( 5 - 10 ) 

|x-yi||x-y 2 | 

It is possible to write (4)h^ using integrals over the matter if this function, G, is used. Eq.(5.10) 



has solutions' 2 ' 1 , 



where 



G(x,yi,y 2 ) = 


- ln(n + r2 ± ri2), 


(5.11) 


n = 


l x -yil> 




r 2 = 


|x-y 2 |, 


(5.12) 


ri2 = 


|yi -Y2|- 





Note that ln(ri +r2 — ru) is not regular on the interval between yi and y 2 , while In(ri + r2 + ri2) 

i« rpcrnlnr nn ihp mat.W Ttm« wp n«p lnffi -I- rn -l-fin"! n« a firppn fnnr+inn TT«incr tViis fnnrtinn 



UUjj and U^Uj are rewritten as 



h l dx^ 1 (y ]x — yi | ^ (,/" |x — y 2 ] 

= / ^,^P( yi )p(y 2 )^( |x _ yi ; x _ y2| ) 
=A/ iaf / (i 2/2p(yi)p(y2)— — — 7 ln(n +r 2 + n 2 ), 

= [ d 3 y 1 d 3 y2p(yi)p(y 2 )-^— 1 ( ] n 1) 

J dy\dy{ M* - yi||x - y 2 |/ 

f d 2 

J ay\ dyi 



(5.13) 



Thus we can express (4)h^ using the integral over the matter as 



iSm " ^ Al ) " 3 (*fe " ^' Al2 ^ ln(ri + r2 + ri2) ' 



(5.14) 



where we introduce 

d 2 



Ai 
A12 



dy\dy\ ' 
9 : 



(5.15) 



QykQyk 



Using relations A/ Iot |x-y| = 2/|x-y| and A/ to |x-y| 3 = 12|x-y|, {A)h^\ (4)^ and (4)4^ 
are solved as 

(4)^ = 2^ j(pv^-y\d 3 y + 2^ J J (pv^ - y\d 3 y + ^ J p|x-y|A, (5.16) 
= ^ /U + 3P-a,x-y,3 y _ V ! ^: P :f' 2 ) ,y. (5,8) 



In total, we obtain 



(4) ft ii - (4)%- 




(5.19) 



5.2. Treatment as a boundary value problem 



The above expression for (4) /if/ is quite interesting because it only consists of integrals over 



the matter. However, in actual numerical simulations, it will take a very long time to perform 
the direct integration. Therefore, we also propose other strategies where Eq.(2.57) is solved as 
the boundary value problem. Here, we would like to emphasize that the boundary condition 
should be imposed at r(— |x|) |yi|, \y2\, but r does not have to be greater than A, where A 
is a typical wave length of gravitational waves. We only need to impose r > i?(a typical size of 
matter). This means that we do not need a large amount of grid numbers compared with the 
case of fully general relativistic simulations, in which we require r > A ^> R. 

First of all, we consider the equation 



Since its source term behaves as 0(r 6 ) at r — > 00, this equation can be accurately solved under 
the boundary condition at r > R as 




(5.20) 




(5.21) 



where 




(5.22) 



Next, we consider the equations for (4)^ , (4)h\j and (4)^ • Noting the identity, 



o = -<WY, = (mfiA .■„■ + Afi^P - (oUi) 4 



(5.23) 



we find the following relations; 



|p|x-y|d 3 



fyl—Lwy, 



p|x-y| 3 d 3 y = 3 / d a y 



l x -y| 

puv ^'^^i"^ + f 4P + ^ 2 - - ^) ) l x - yl 



(5.24) 



Using Eqs.(5.24), (4)h^\ (4)^y^ and (4)h\J^ in Eqs.(5. 16-18) can be rewritten as 



(4) 



h 



(C) 



2 / (pvi) 



x - y .3 



l x ~ y 



d 6 y + 2 / (pu*) 



l x -y| 



M .x k -y k « 



d 6 y, (5.25) 



^ = Ad^i j pv v — ^ — y + y ^ } v^\ d y 

1 \ d f - ^) ,3 9 f r>/ (x i - y*) o 

+ o i 7TT / P 1 ^ + 7T7 / P 1 ~T d V 

2 yox 1 J |x — y| oxJ J |x — y] 

ox J |x — y ox-> J x — y 



(5.26) 



and 



( 4 )/l 



(V) 



pv 2 + 3P 



_ pC/\ ^ - //•' ; 



2 / |x-y| 



d 6 y 



d 

dxi 



pv 1 + 3P 



d y 



2 J x-y 



2 r(pv 2 + ?,P- P U/2) 

" 3^ y ^ ~ 



<Py, 



(5.27) 

where P' = P + pv 2 /4 + pU^y 1 /A. From the above relations, ^h^\ (4)h\f and (4)^-J^ become 



(4) ^f = 2(^(3)^' + ^ (3) P* - Qy) + - ^ (3) P fc ) , 

w% - ± dx i dx j{ Vki xx ZVk x +v ) + 3 di A x & n 2 J 

+ ^H*" - + - vr »)}, 



where 



A flat (3) Pi = ~^pv\ 

AflatQij = -47r[^'(p^)- +X i (pv j )-y 
&flatQ (I) = ~^[pv 2 + 3P- ± P U), 

^flatQ { P = ~^(pv 2 + 3P- l - P u). ' 



w l v J , 



A/fc^ = -4^ 

AflatV^ = -47T/CWVV, 

A //at y^) = -47rp(uV) 2 
A /to \/( p ) = -47TP', 
A //at y/ P) = -4ttPV, 
A/^ (pf/) = -4^, 
A //at ^ f/) = -AnpUixi. 



(5.29) 



nan De derived from the above potentials which satisfy the 



Therefore, (4)hy\ (4)^ and (4)^)^ can 
Poisson equations with compact sources. 

We note that instead of the above procedure, we may solve the Poisson equation for ^hij 
carefully imposing the boundary condition for r ^> R as 



3 k 

-n K 

4 



_ ZtjVtjM 2 ) 4- -n i n :j n k n l n m 4- —A. 7"^ 2 ) _ -S- n k n l n m T^ \ 
g n n n 1 kll + ^nn n nn i klm -+- ^o tJ n i m g o y n n n -* feZm j> 



(5.30) 



+ I 3 n kj + jki ^ ~~ s^Sjkk + nJ Sikk) 

+ 2n k n\n l S jM + n 3 S iU ) + 2n l n 3 n k S k u + -^n^jJ + 0(r~ 3 ). 

It is verified that 0(r _1 ) and 0(r~ 2 ) parts satisfy the traceless and divergence-free conditions 
respectively. It should be noticed that (4)/iy obtained in this way becomes meaningless at the 
far zone because Eq.(2.57), from which (4)/iy is derived, is valid only in the near zone. 



6. Conserved quantities 



The conserved quantities are gauge- invariant so that, in general relativity, they play impor- 
tant roles because we are able to compare various systems described in different gauge conditions 
using them. From the practical view, these are also useful for checking the numerical accuracy 
in simulations. Thus, in this section, we show several conserved quantities in the 2PN approxi- 
mation. 

6.1. Conserved Mass And Energy 

In general relativity, the volume integral of the mass density p does not conserve, and instead 
we have the following conserved mass; 

M * = J P*d Sx - (6- 1 ) 
In the PN approximation, p* is expanded as 



P* = P 



+ 317) 

+ + \v 2 U + ^U 2 + 6 (4) V + ( 3 )/V) + (6)5* + 0(c- 7 ) 



(6.2) 



where ( 6 )<5* denotes the 3PN contribution to p*. 

Then, we consider the ADM mass which is also the conserved quantity. Since the asymptotic 



behavior of the conformal factor becomes 



i> = i + 



Madm 
2r 

the ADM mass in the PN approximation becomes 



+ 0(1), (6.3) 



M ADM = ~ J &fiati>d 3 

= J d 3 xp[{l + {v 2 + e + lu) + (V + ^v 2 U + v 2 e + -J + \ue + \u 2 (6.4) 
+ 5 (4) ^ + 2 {3) PiV*) } + ^ ( (3) iy (3) iy - l(3)K 2 ) + {6) 5adm + 0(cT 7 )] , 



where ^Sadm denotes the 3PN contribution. 



Using these two conserved quantities, we can define the conserved energy as follows; 



/ -' =M ADM - M* 

+ + 3v 2 U + v 2 e + -v 2 + hje - hj 2 - (4) V + (3 )/W) } 



+ 



1 



f 67rp 

iEn + Eipn + E2PN + 



(3) A ij{3) A ij - 3(3)^) + ((6)^M - (6) 5 *) + °( c ? )] 



(6.5) 



We should notice that the following equation holds 



J (3)4(3)4'^ = -8t y pv\ 3) (3id 3 x + J (l(3)K 2 + 2f/ (3) ir)rf 3 a;, (6.6) 



where we use the identities derived from Eqs.(2.50) and (2.51) 



{s) Pij i 3)Pijd 3 x = -ltor J pv\ 3) (3 i d 3 x + J( {3) K 2 + 2U {3) K-3U 2 )d 3 x, 
J {3) (3 hj{3) (3 j , l d s x = j{ {i) K 2 + W {i) K + W 2 )d 3 x. 

Using these relations, we obtain the Newtonian and the first PN energies as 



En = J p(\ v2 + £ -l U ) d3 ^ 



(6.7) 



(6. 



and 



E\pn = I d x 



5 4 5 2 



E\pn can be rewritten immediately in the following form used by Chandrasekhar [261 ; 



E 1PN = J d 3 xp + ^v 2 U + v \ £ + ^) + 



2Ue - -W 



(6.10) 



where qi is the first PN shift vector in the standard PN gauge(i.e., = 0) and satisfies 



E2PN is calculated from the 3PN quantities ( 6 )£* and ^Sadm- (6)^* becomes 

(6)^ + ^v'U + , 2 (5 (4) ^ + | f/ 2 + | (3) /W - X) 

+ 6 (6) ^ + 15Z7 (4) ^ + ^t/ 3 + 7 (3) /VC/ + ^hijvV + ^ (3) A(3)A + (5)/V, 

and we obtain 
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(6)Py -P 



+ (4 )MV' + 2 (3 )A^ i {2^ 2 + £ + ^ + y f/ } + 2 (5)A^ + (3)A(3)A 

Making use of relations (tythjj = and (6)^ij,j = 0, we obtain 

(6)Madm = J d 3 x^p^p 

+ ~Ll ^vfaAlsA - !(3)*(5)*) ( 6 - 16 ) 

+ 3^/ ^((3)4(3)4- -| 3) ^ 2 ), 

where we assume (6)^ij — ► 0(r _1 ) asr — > oo. Although this assumption must be verified 
by performing the 3PN expansions which have not been done here, it seems reasonable in the 
asvmntoticallv flat SDacetime. From i^Madm and ^M*. we obtain the conserved enerev at 



(6.12) 



( 6 )M* = J p^S+dPx. (6.13) 
The Hamiltonian constraint at 0(c~ 8 ) becomes 

^flat{&)i> - {4)hij(4)i>,ij - -^{<S)hijU,ij 

= -7^(2(4)hkl,m(4)hkm,l + (4)hkl,m(4)hkl,m) (6.14) 

1 / ~ ~ 2 \ 1 / ~ ~ 2 2 \ 

- ^(6)Pi> - l{(3) A ij(5) A ij ~ ^(3) K (5) K ) ~ {(3) A ij(S) A ij ~ 3 (3) K ), 

where we define (6)Py as 

; 6 + t; 4 (e + ^ + yC/) + ^ 2 {f + ^) + 9 (4) V ~ 2X + 20C/ 2 } 
+ e(5(4)^ + ^ 2 ) + 5( 6 )^ + 10^(4)^ + ^ 3 (6-15) 



the 2PN order 



E2PN =(6)MadM - (6)-^* 

= [d 3 xp\ 1 -^v G + v 4 (e + - + ^u) 
J L 16 V p 8 / 

+ - 2 {4(4)^ - X + H u(e + + ft/ 2 + 
+ e (5 (4) V + - ( 6 )^ " 5t/ ( 4) V - 

+ l^hijvV + 2 {s) p i v i (e + + 3C/) + (5) /V + ^(3)A(s)A} 



+ 
+ 



1 
1 

327T 



d 3 y uL ) A ij{3) A ij -- {3) K 2 



When we use the relation, J d^xp^ip = — J (fixlIA^ip, we obtain 

E 2PN = J S>xp [IV + , 4 (e + - p + f £/) 

+ t; 2 {4 (4) ^ - X + 617 (e + ^) + yf/ 2 + | (3 )/W} 
+ e(5 (4) ^ + ^ 2 )-yf/ (4) ^-^ 3 

5C/} +(5)A^ + ^(3)A(3)A 



+ -(4)MV+ 2 (3)A^{( 



+ 



8tt 



(6.17) 



(6.18) 



6.2. Conserved linear momentum 

When we use the center of mass system as usual, the linear momentum of the system should 
vanish. However, it may arise from numerical errors in numerical calculations. Since it is useful 
for investigation of the numerical accuracy, we mention the linear momentum derived from 



P, 



, - - lim <b[Kijn j - KrAdS 



1 

8^ 

:— lim d>[ ^A ijn J - -Kn l )dS, 
87T r^oo JV l 3 3 / ' 



(6.19) 



tirh oro \- h a on t* f • 1 r*c m+ orrvo 1 o o t*o +■ o r^^i ror o cnliara r\-f r»/~\n o^on-f - <r* Qm r»o f Iid o rim t~i^ V\ on qi n f^"v 



of Aij is determined by 



1 2 \ 

(3)^j = 2 ((3)A,i + (3)/5j,< - g^-(3)A,iJ + 0(r- 3 ), 



and 



(5)A? = ^ ((5)Aj + (5)/?j,t - j^y(5)A,j) + (5)^5 T + 0(r 3 ), 
the leading term of the shift vector is necessary. Using the asymptotic behavior 



(3)A 



7Zv InVI, 



2 r 2 r 



+ 0(r- 2 ), 



the following relation is obtained 



2 \ . 



where k = J pv i d^x. Therefore the Newtonian linear momentum is 



P N % = l d 3 xpv i 



Similarly the first PN linear momentum is obtained as follows; 



Pipn 1 = / d 3 xp 



i f „,2 



P 



r' (('" + £ + 6U + —J + (3) A 



P2PN 1 is obtained by the similar procedure. 



7. Summary 



In this paper, we have developed the PN approximation in the (3+1) formalism of general 
relativity. In this formalism, it is clarified what kind of gauge condition is suitable for each 
problem such as how to extract the waveforms of gravitational waves and how to describe 
equilibrium configurations. It was found that the combination of the conformal slice and the 
transverse gauge is useful to separate the wave part and the non-wave part in the metric variables 
such as hij and if). We also found that, in order to describe the equilibrium configuration, the 
conformal slice is not useful and instead we had better use the maximal slice. Although we 
restricted ourselves within some gauge conditions in this paper, we may try to use any gauge 
condition and investigate its property relatively easily in the (3+1) formalism, compared with 
in the standard PN approximation performed so far 1241 . 

We have also developed a formalism for the hydrodynamic equation accurate up to 2.5PN 
order. For the sake of an actual numerical simulation, we carefully consider methods to solve 
the various metric quantities, especially, the 2PN tensor potential (4)/iy- We found it possible to 
solve them by using standard numerical methods. Thus, the formalism developed in this paper 
will be useful also in numerical calculations. 

Here, we would like to emphasize that from the 2PN order, the tensor part of the 3- 
metric, 7^, cannot be neglected even if we ignore gravitational waves. Recently, Wilson and 

[27] 

Mathews presented numerical equilibrium configurations of binary neutron stars using a semi- 
relativistic approximation, in which they assume the spatially conformal flat metric as the spa- 
tial 3-metric, i.e., jij = 5{j. Thus, in their method, a 2PN term, hij, was completely neglected. 
This means that their results unavoidably have an error of the 2PN order which will become 
~ (M/R) 2 ~ 1 — 10%. If we hope to obtain a general relativistic eqiulibrium configuration of 
binary neutron stars with a better accuracy (say less than 1%), we should take into account the 
tensor part of the 3-metric. 

In section 3, we used several slice conditions and investigated their properties, but, as for 
the spatial gauge condition, we fix it to the transverse gauge for the sake of convenience. It is 
not clear, however, whether this is the best gauge condition in numerical relativity. In numerical 

t*o1 o 4" i~k n4~^ r 4" It c> c Vi 1 ff -\ ropf rw Til en re o s. rori r 1 m t~i /"it* 4* o n+ rr\] o 4" r\ "vorl 1 1 r»o 4" V> a /~>r\r\v*r\ l n a 4" q oIidov uzo 4*o 1 1 



to choose the appropriate condition, the coordinate shear in the spatial metric will continue to 
grow, and as a result, the simulation will break down. The minimal distortion gauge which was 

[28] 

proposed by Smarr and York is a candidate which may efficiently reduce the coordinate shear. 
Even if we use this gauge condition in the PN analysis, equations for u\hij, f$\hij, (3) A, (5)$ 
and (6) A remain unchanged, but higher order terms of hij and may slightly change. If we 
investigate the effects due to the difference, we may be able to give some important suggestions 
about the gauge condition appropriate for numerical relativity. 

In this paper, we mainly payed attention to the PN hydrodynamic equation in the (3+1) 
formalism to describe the final hydrodynamic merging phase of coalescing binary neutron stars. 
This formalism, however, may be useful also for gravitational wave generation problems in 
the inspiraling phase, which is usually performed using the harmonic gauge' 4 ' 51 . Since we can 
choose any gauge conditions in the (3+1) formalism, it may be possible to reduce some lengthy 
calculations which appear in previous works using the harmonic gauge' 4 ' 51 . This problem is 
interesting and important as a future work. 
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APPENDIX A 

Calculation of ^hij 

We make use of the transverse property of r^, which is guaranteed by the transverse gauge 
condition, in order to obtain (5)/iy. Using the following identity 

fA\Tin = (lA\TSUXp\ u . (A.l) 



Eq.(2.58) can be rewritten in the surface integral form 



1 d f 

(5) h ij = -^qI J (4)T ik x J n k dS. (A.2) 

Thus, we only need to estimate terms of 0(r~ 3 ) in (4)%-, which come only from the shift vector 
in the conformal slice as 

(3)& = ^2 ( nj2 ij + \ nji ij + \ nii kk - \n i -n?n k i jk + ^-J {3) Kd 3 x^j + 0(r- 3 ), (A.3) 



where 



Zij(t) = -A J pv i y j d 3 x. (A.4) 
Here, note the following relations as 



dt 
and 



9 Zij = -2/y, (A.5) 



{3) Kd 3 x = 2nl kk . (A.6) 



Therefore, the relevant terms of (4)Ty for the surface integral become 

2 



^ [{ -3 Iij + 3 (/^n V + Ij k n k rf) - - 4 fc raV + y 4znW™ V } (A.7) 



Thus we obtain fry at the 2.5PN order, 



hij = J ((3)A,jfc + (3)4,i - I^Ai)^"*^ 



5 V [ ) - 



(A. 



This derivation seems fairly simple owing to the gauge condition. Thus, it is expected that 
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APPENDIX B 



In this appendix, we briefly comment on a method to solve the Poisson equations for (fyfii 
and (5)A, i.e., Eqs.(2.51) and (2.53). Since the source terms of them have terms such as —U,i 
and —2(4)^ which behaves as 0(r~ 2 ) at r — > oo^, it seems that a technical problem arises in 
solving these equations in numerical calculation. However, this is easily overcome in a simple 
manner. We consider the case of the maximal slice for simplicity, but other cases may be treated 
similarly. 

First of all, we write ( 3 )/3j and ( 5 )/3j as, 



(3) A = -4(3)-P< + 7}X,i, 



(b.i; 



(5)A = -4(5)-P» + 9(4)X,t, 



where x an d (3)-Pi satisfy Eq.(4.4) and the first equation of Eqs.(5.29), respectively. (5) Pi and 
(4)X satisfy the following Poisson equations; 



P\ 



&flat(5)Pi = -4vrpp^ 2 + 2U + e + -) + (3) a] + 2^(3)^ - ^(UU) 



Aflat(4)X = -4(4)^- 



(4)X 



can be written as 



(4)X 



-/ 



p 4 |x-y|rf 3 y, 



^((3)A^,/),i , 

(B.2) 
(B.3) 



where 



P4 = p(^+£+^). 

From Eqs.(4.4) and (B.3), x,i an d (4)X,i become 



X, 



= - d yp- 



i x ~ y 



(4)X,» = 



^ 3 1/P4 | x _ y | = ~2x l (4)^ + ( 4 )7/i, 



(B.4) 



(B.5) 



f Note that in the Newtonian limit, —U,i is 0(r 4 ), but at the 1PN order, it becomes 0(r 2 ) because 



where rji and ^rji satisfy 

AflatVi = ~^px\ 



(B.6) 

Aflat(4)Vi = -4:7rp 4 x\ 



Hence, 

(3)A = -4 (3 )Pi -l(x l if -m), 

\ V 7 (B.7) 

(5)A = -4(5)-P» - g ( 2a;, (4)^ - (4)^)- 

Since the source terms of the Poisson equations for (3) -Pi, (5)Pi> ?7i an d (4)^ behaves as 0(r _n ), 
where n > 5, at r — > 00*, these vector potentials can be accurately obtained by solving the 
Poisson equations for them under appropriate boundary conditions. Thus, there is no difficulty 
to obtain (3)/% and (5) A- 

Finally, we note that the above method is not unique prescription. For example, (3)/% in the 
first PN approximation"'" may be expressed as 

(3)A = -4 (3 )Pi + \{{ { z)P k x k ),i - X2,i}» (B-8) 

where %2 satisfies 

A flatX2 = -^pv i x i . (B.9) 



* Note that the non-compact sources of the Poisson equation for (5)Pi may be regarded as 0(r 5 ) in the 
2PN approximation because U is 0(r~ 3 ) in the Newtonian order. 
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